Adaptation of Autocatalytic Fluctuations to Diffusive Noise 
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Evolution of a system of diffusing and proliferating mortal reactants is analyzed in the presence 
of randomly moving catalysts. While the continuum description of the problem predicts reactant 
extinction as the average growth rate becomes negative, growth rate fluctuations induced by the 
discrete nature of the agents are shown to allow for an active phase, where reactants proliferate as 
their spatial configuration adapts to the fluctuations of the catalysts density. The model is explored 
by employing field theoretical techniques, numerical simulations and strong coupling analysis. For 
d < 2, the system is shown to exhibits an active phase at any growth rate, while for d > 2 a kinetic 
phase transition is predicted. The applicability of this model as a prototype for a host of phenomena 
which exhibit self organization is discussed. 
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I. INTRODUCTION 

There is a growing interest in the dynamics of cat- 
alytic systems with diffusing reactants |IJ . These models 
has been considered in the theory of population biology 
[2 1, chemical kinetics || and physics of contact process 
\i\ as well as magnetic systems ||. In the most simple 
case, where agents undergo only birth (autocatalytic re- 
production) and death (spontaneous annihilation), the 
total growth rate, i.e., the difference between the typical 
rates for these two processes, is the critical parameter for 
the system. While negative growth rate implies exponen- 
tial decrease in the number of particles toward extinction 
of the "colony" , positive rate gives exponential growth. 
Usually, the number of reactants saturates to some con- 
stant value which reflects the finite holding capacity of 
the environment. The most simplified mathematical de- 
scription of this process is given by the continuum Fisher 
equation ||: 



db(x,t) 
dt 



DV 2 b(x, t) + ab{x, t) - Xb 2 (x, t) 



(1) 



where b(x, t) stands for the density of reactants, a is the 
total growth rate, and —Xb 2 is introduced phenomeno- 
logically as the "minimal" nonlinear term which leads to 
saturation at positive growth rate. Since the density b 
is positive semi-definite, at negative growth rate (a < 0) 
there is only one steady state, the absorbing state, where 
6(x, t) = everywhere. At positive tr, this state be- 
comes unstable, and the system flows into the uniform 
state b — a/X. In this simplified framework the diffusion 
is irrelevant to the steady state, and only governs the 
dynamical approach to it, an effect which has been con- 
sidered at Q . It turns out that the typical invasion of the 
unstable phase by the s table one is in the form of Fisher 
fronts (of width w ~ y/D/tr) which propagate with ve- 
locity v > 2\jDo. At the stable state, any small fluc- 
tuation with wavelength k decays as exp[— (er + Dk 2 )i\. 
The phase transition from the inactive to the active state 
takes place at a = 0. 



Equation (|l]) describes the continuum limit of many 
underlying discrete processes. A typical example is a 
system with particles diffusing (random walk), annihi- 
lating (B + B — > 0) and reproducing autocatalytically 
B — > 2B. The discrete nature of individual reactants 
and their stochastic motion introduce a (multiplicative) 
noise, which may dominate the evolution of the system 
and violate the predictions of (|l|). For the above men- 
tioned and similar processes, it turns out that, at low 
dimensionality (d < 2), the extinction phase is stable 
even at small, positive growth rate, and the transition 
from active to inactive state falls into the equivalence 
class of directed percolation || (Reggeon field theory). 
Moreover, it has been conjectured |t]] that any transition 
with single absorbing state falls into the same equivalence 
class, unless some special symmetry or conservation laws 
are introduced, as the even offspring case considered by 
Grassberger ||, and Cardy and Tauber ||. 

Recently, a new type of active-inactive phase transition 
has been introduced, for the process (1) in the presence of 
quenched disorder in the relevant term, a{x). It has been 
shown by Janssen Jl0[ , that the renormalization group 
(RG) flow of that process has only a runaway solutions 
in the physical domain (due to the "ladder" diagrams 
which changes the effective mass of the free propagator), 
hence the phase transition is not of directed percolation 
type. Nelson and Shnerb jLl| , using the continuum ap- 
proximation, showed that the local growth is related only 
to <r(x) at the vicinity of the domain, i.e., localization of 
Anderson type Jljj takes place and the extinction transi- 
tion is given by the effective growth rate of small, local- 
ized islands. Although the effect of intrinsic noise due to 
discreteness fluctuations has not been considered in [ ll) , 
it seems reasonable that the actual transition takes place 
when the time scale for tunneling between two positive 
growth islands is smaller than the time scale for absorb- 
ing state decay of a single "oasis" . 

If the disorder is uncorrected, cr(x, t) with 
(ct(x, tWx't')) = A<5(x - x')<5(t - t'), the linear part 
of Eq. (pi) also governs the statistics of a directed poly- 
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mer on a heterogeneous substance, where the time in 
(jlh is identified with the polymer's preferred direction 
Jl3[ . With the Cole-Hopf transformation, this problem 
is mapped to the noisy Burger's process [Q (KPZ sur- 
face growth fill). In contrast with the "localization" in 
the case of static disorder, uncorrelated environment in- 
duces superdiffusion of the reactants, where the "center 
of mass" of the population wanders in space as r oc fi, 
with ( > 0.5. The effect of intrinsic stochasticity due 
to discretization is, again, limited, and the statistical 
properties of the eigenenergies of the directed polymer 
problem determines the extinction transition. This is 
also the case when the system under quenched disorder 
is subject to strong convection, as has been shown by 

In this paper, we consider the case of diffusive disorder. 



B + A 



2B + A 



when both B and A undergo diffusion with rates and 
D a , respectively. The mortal agent, B, dies at rate /i, 
and proliferates in the presence of the (eternal) catalyst 
A. The continuum description for this process is given by 
the "mean-field" (rate) equations for the densities a(x, t) 
and b(x, t): 



db(x,t) 
dt 



= D b V 2 b(x,t) - fib + Xab 



da(x, t) 
dt 



D a V 2 a(x,t). 



(2) 



As t — > oo, a flows into its average a, thus the effective 
mortality rate for b is given by m = /i — Xa (the mortal- 
ity rate turns out to be the "mass" of the effective field 
theory, hence denoted by m). For positive "mass" the 
b population decays exponentially while negative mass 
implies exponential growth. The active-inactive phase 
transition takes place at m = 0. One observes that equa- 
tion ^) is obtained from ([!]) by dropping the non-linear 
term and replacing a by Xa — u, accordingly, at long 
times the process introduced in (J2J) is the linearized form 
of (Q) with the proliferation rate, a, fluctuates "diffu- 
sively" around its mean Xa — fi. As for this system the 
disorder is not static but is correlated, it somehow in- 
terpolates between the above mentioned models and one 
may wonder weather it leads to localization of the re- 
actants or to superdiffusion flqj . It turns out that the 
reactants may adapt themselves to the environment and 
the colonies are localized on the diffusing islands. More- 
over, these correlated fluctuations due to the stochastic 
motion of individual reactants will change the character 
of the transition; the transition point is pushed to neg- 
ative values of m and the b reactants survive below the 
"classical" threshold. Some of our results, along with 
a numerical study of the transition, are summarized in 
previous publication Jl7[. 



II. STRONG COUPLING ANALYSIS 

The basic intuition beneath the phenomenon we de- 
scribe is in the concept of adaptive fluctuations. Let us 
take a look, first, at the case of frozen As, where we have 
random, quenched, growth rate as in [ ]TT[ . The linearized 
continuum equation then takes the form 



db 
dt 



DV 2 b-m(x)b 



where m(x) = jj,— Xua{x), and ua{x) is the random con- 
centration of the catalyst A. The system is in its active 
phase if the linearized evolution operator 

L = DV 2 - m(x) 

admits at least one positive eigenvalue, and its localiza- 
tion properties are almost determined by the correspond- 
ing eigenfunctions. Since the same operator governs the 
physics of quantum particles in random potential, one 
may use the known results |l8| of this field, i.e., that 
in low dimensionality or strong disorder all the wave- 
functions are exponentially localized and the diffusion 
becomes irrelevant on large length scales 19 1. Accord- 



ingly, the system may be in its active phase at localized 
islands even if the average m is positive. 

In our case, however, the catalysts diffuse, and these 
colonies survives only if the reactants cluster is able to 
"trace" a specific catalyst or to find some other wander- 
ing island. This implies the significance of the system 
dimensionality: while two typical random walkers (such 
as the catalyst and the reactant) encounter each other in 
finite time for d < 2, they will (typically) never collide 
for d > 2. One may expect, accordingly, that below 2d 
quantization induced fluctuations are much more domi- 
nant than above two dimensions. 

Consider one "frozen" catalyst at the origin. The ef- 
fective growth rate in the vicinity of the origin is positive, 
i.e., m(r < R) — rrii n < in a region of typical catalyst 
size R around it. In the "desert", out of this island, the 
"mass", m out , is positive. The colony is then localized 
at r — 0, with growth rate \mj n \ and an exponentially 
decreasing tail into the desert In the continuum 

approximation, the time dependent profile of the tail is 
given by p0|: 



b(r) 



,\r 



t—ryr, 



t/D 



(3) 



and the tail front, which is the size of the reactants 
colony, moves away from the origin with typical veloc- 
ity v ~ m in yj D /m ou t. If the catalyst is moving, the 
colony will die only when the A molecule detaches from 
the B colony. This, however is almost impossible for a 
diffusively moving A since the colony's front moves bal- 
listically. As this argument involves only one catalyst it 
is independent of system dimensionality Thus, at strong 



2 



Log B Concentralion 




FIG. 1. The profile of a B island as a function of time as it 
follows the random motion of an A agent. The cross-section 
of the island is taken through the current location of the A 
agent. The inset shows the time evolution of the logarithm 
of the height of the B concentration at the point at which 
A is currently located (solid line). The B colony is seen to 
grow, although the average growth rate over the entire space 
is negative (ua is extremely low since there is only one in the 
whole simulation space, thus Xua — M ~ — n)- The dashed 
line shows the exponential growth according to Eq. (^). 

coupling some localized islands are active in any dimen- 
sion, in contradiction with the mean field prediction of 
extinction at positive average mass. This mechanism is 
illustrated in figure ([!]), which manifests the ability of 
the B colony to adapt to the location of the moving cat- 
alyst. We stress that this adaptive skill is solely due to 
the "dumb" diffusion and multiplication of the reactants, 
and thus is an emergent self-organized feature. 

There is, however, a possibility for a different scenario, 
the weak coupling limit, where the local properties of the 
system do not support the formation of colonies in the in- 
active phase. For a B molecule having a spatial overlap 
with an A catalyst, the multiplication time is propor- 
tional to A. If this time is much larger than the relevant 
hopping rate, the typical birth event is singular and no 
colony is formed. If, furthermore, the decay time for a 
B particle in the absence of A is much smaller than the 
typical time to find a new catalysts one may expect that 
the system is in its inactive phase, unless some global, 
collective effect turns this local analysis void. 

Before looking for global effects, let us try to consider 
the strong coupling limit more carefully. Consider a sin- 
gle A agent located at the point Tq, and the "island" of 
B reactants which surrounds it. Keeping the A station- 
ary and working on a d-dimensional lattice (with lattice 
constant / , the growth rate is jj and the hopping rate is 
the following equation holds for the concentration 
of S's: 



n B 



(4) 



where we assume D B l d ~ 2 -C A (thus the very steep 
slope suppresses the effect of diffusion returning inwards). 
Consequently, 



dlog[n B (r ,t)} (\-2dD B 



Of 



I 2 



(5) 



Now consider a hopping event of the A, by a single lat- 
tice spacing. Measuring n B (r (t) , i), at the new A site, 
it reduces by a factor of 



(6) 



The rate at which the hopping events occur is 2 S A , ac- 
cordingly the rate equation (||) is modified to: 



dlog [n B (r jt),t)} 

at 

A - 2dD B 
I 2 



2dD A f, d _ 2 D B 

m — — log I r - 



I 2 



A ) ' 



(7) 



where we have assumed that there is enough time be- 
tween hopping events so that the island shape stabilizes 
to the long time behavior (Q), namely DaI^ 2 <C A. 
Equation (^) has the solution: 



»B 



(\-1dD B 2dD A f,d-20 b \\ 

(ro(t),*)~e(— i 3 "~ i^M' ^)) t . ( 8 ) 



This shows that in the strong coupling regime the ex- 
ponent is positive and the number of reactants grows 
exponentially, independent of the dimension and of the 
catalyst density. The inset of Fig. (Q) shows the fit of 
the expression (||) to the numerical results on a lattice. 



III. WEAK COUPLING 

In order to consider global effects of spatial fluctua- 
tions, let us write the Master equation for the probability 
Pnm to find m reactants and n catalysts at a single point 
(with no diffusion) 



dP n 



dt 



= - H [mP nm - (m + l)P„, m+ i] 

- X[mnP nm -n(m- l)P n , m -i]- (9) 



Following plj we define a set of creation-annihilation 
operators, 

o + |n, m) = \n + 1, m) b + \n, m) = |n, m + 1) 
a |n, m) = n \n — 1, m) b \n, m) — m \n, m — 1) 

and a wavefunction 

* = P n ,m \n, m) . 

The Master equation then takes the Hamiltonian form: 



dt 



-Hi- 



wit h 



3 



<e-i> 



+ iiHk - h] + -[a+a.b+h - a+ ai b+b+b 



(10) 



Where i runs over all lattice points, and the sum, (i — e), 
is over nearest neighbors. 

Shifting the creation operators to their vacuum expec- 
tation value a + — > a + 1 and b + — > 6 + 1 , the value of the 
catalyst density to its average, p-a — > a + n a , and finally, 
-pb — > 6 the evolution operator takes the following form 
in the continuum limit: 

H = J d d x [ -D b bV 2 b - D a aV 2 a + isbb 

- \bb(a+l)(b+l)(a + n a )\ (11) 
The action is simply: 
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d x a—a + b—b + H 



at 



at 



(12) 



Note that the coefficient of 66, which plays the role of 
mass, is given by m = ji — Xn a . 

Now this system may be analyzed using the standard 
renormalization group (RG) technique. We impose a 
change of scale 



x ■ 
t - 



sx 
-> s z t 



(13) 



A — * A/s 



where s is the renormalization group scale factor. The 
renormalization flows of the parameters of the action (|l^ ) 
are given by their naive dimensionality and the correc- 
tions from the diagrams shown in Fig. (Q), using the 
basic vertices as in Fig. (^). The flow equations for the 
mass and the coupling constant are given by 



dX 



d ln( S ; 



= e\ + 



A 2 



A" 



2ttD 1 



dm A 2 no 

= 2m — 

2vrL> 



d ln(s) 



and the flow lines are shown in Figs. and (||) for 
d < 2 and d > 2 , respectively. Below two dimensions, 
the Gaussian fixed point is always unstable and the sys- 
tem flows to the strong coupling limit, where adaptive 
B colonies grow indefinitely, as indicated by the negative 
values of the effective mass. At higher dimensionality, on 
the other hand, there is a finite region in the parameter 
space where the trivial fixed point is stable and A flows 
into zero, while higher values of initial coupling constant 
flow to infinity. For a system of finite size L d , the flows 
should be truncated at s = j, and the phase is deter- 
mined by the end point of the flow lines at this s. For 
finite A, equations (|l4| ) take the form: 



d") 



d ln(s) 



dM 
d ln(a) 



£7- 



2M 



7 



1 + M 



Ct"f 

1 + M' 



AA 



M 



where the dimensionless quantities are 7 — — err- 
and a — 2 ^" a . The flow lines for d — 2 are shown in 



figure @ and exhibit a transition due to the finite lattice 
spacing. 



'b - b 
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FIG. 2. Elements of perturbative expansion: the vertices 
corresponding to the action ( p^[ ) 



DA 2 



dm 
d Ms) 



= 2m 



A 2 r_ 
2nD 1 



A 



d-2 



(14) 



DA 2 



where e = 2 — d, D = 



Da+D b 



and A(s) is the upper mo- 



mentum cutoff. Note that, as indicated by naive dimen- 
sional analysis, the Gaussian fixed point {A = 0, m = 0} 
is stable at d > 2, hence there is no perturbative correc- 
tions to rj and z in this regime. If the momentum cutoff 
is much larger than any other quantity of the problem 
one has, 



r/A 



d ln(s) 



A e + 



A 



2ttD 



FIG. 3. Most UV divergent diagrams contributing to the 
renormalization group equations. 
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FIG. 4. The RG flow lines at the continuum limit for 
d < 2. While at short times m grows, the system flows into 
its active phase m < on large time scales. 
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FIG. 5. Renormalization flows for d > 2. The shaded re- 
gion flows into to active phase while at the unshaded region 
the system flows to the inactive phase (m — > oo). 
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FIG. 6. Renormalization flows for d = 2 at finite lattice 
spacing. Unlike (^), here there is a region in parameter space 
(unshaded) where the extinction phase is stable. 

IV. DISCUSSION 

Since the classical works of Malthus and Verhulst , 
it has been recognized that most of the processes in liv- 
ing systems are autocatalytic and thus are characterized 
by exponential growth. In fact, the appearance of an au- 
tocatalytic molecule may be considered as the origin of 
life. In this paper, these autocatalytic system are shown 
to admit self organization in the presence of fluctuating 
environment. The exponential amplification of "good" 
fluctuations in the catalysis parameters prevails, in the 
situations discussed above, the globally hostile environ- 
ment and is robust against the random motion of both 



the reactants and the catalysts. Our result may be inter- 
preted as an indication that "life" (in the above sense) is 
resilient and is able to adapt itself to the changing envi- 
ronment. The applicability of this model ranges from bi- 
ological evolution (where the environment is the genome 
space) to the role of enzymes in chemical reactions and 
even in social or financial settings. 

More realistic models, however, should take into ac- 
count the depletion of resources by the catalytic process 
and the finite carrying capacity of the substrate. Al- 
though the model discussed above is relevant at time 
scales which are small in comparison with the mean time 
for consumption or saturation, the stable fixed point of 
the system may be different. In particular, on a uniform, 
inexhaustible, autocatalytic substrate with finite carry- 
ing capacity the discreteness induced fluctuations have 
been shown ]9| to decrease the effective growth rate, and 
to give a directed percolation type transition at d < 2. 
The competition between this effect and the effect of 
adaptive fluctuations will be considered elsewhere. 
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